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ABSTRACT 

In astrophysical magnetohydrodynamics (MHD) and electrodynamics simulations, numeri- 
cally enforcing the V • -B = constraint on the magnetic field has been difficult. We observe 
that for point-based discretization, as used in finite-difference type and pseudo-spectral meth- 
ods, the V -B = constraint can be satisfied entirely by a choice of interpolation used to 
define the derivatives of B. As an example we demonstrate a new class of finite-difference 
type derivative operators on a regular grid which has the V • -B = property. This principle 
clarifies the nature of V B ^ errors. The principles and techniques demonstrated in this 
paper are particularly useful for the magnetic field, but can be applied to any vector field. 
This paper serves as a brief introduction to the method and demonstrates an implementation 
showing convergence. 
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1 INTRODUCTION 

As originally laid out by |Brackbill & BamesH1980) , failing to obey 
the V • i? = constraint in magnetohydrodynamics may lead to 
numerical instability and unphysical results. This issue has been an 
issue which has attracted much attention in computational astro- 



physics (ex: Brackbill & Barnes 


1980||BaIsara & Kim|2004||Price| 


|2010| |Dolag & Stasyszyn||2009 


1. To elucidate what V • B = 



means, specifying the manner in which B is represented is essen- 
tial. In a numerical method, the vector fields are represented by a 
discrete set of values. Two classes of discretizations are popular 
in astrophysical applications, finite-volume and point values. Fi- 
nite volume discretizations store the volume-average of the field 
over some cell. These volume averages constrain the possible di- 
vergence of a vector field interpolating these values, and hence the 
Constrained Transport method l |Evans & Hawley|1988| > can be ap- 
plied to conserve this divergence throughout the simulation. How- 
ever, when the magnetic field is represented by point values, the 
divergence of the interpolated field is not constrained by the point 
values, so some extra freedom exists. Two classes of approaches 
have been used. The first class is to admit V-B 7^ errors, and then 
attempt to manage the consequences. Methods of this type include 
the 8-wave scheme (Powell" 1994'; 'Powell et al.'1999), and diffu- 
sion method ( Dedner et al. 2002 1. The Smoothed Particle Hydrody- 



namics schemes of 



IPrice & Monag haiil j2004a|b| |2005t and |B0rve,| 
|Omang & Trulsen] 2001^ ; |Dolag & Stasyszyn "l |2009 ^ also fall into 



this class, as the former uses a formulation of the MHD equations 
which is consistent even in V • B 7^ 0, and the latter removes 



the V -B 7^ contributions to the momentum equation. The sec- 
ond class of methods constrain the derivatives of the interpolated 
field. The projection method, used in finite-difference jBrackbilT] 
|& Bames|1980[ l, and pseudo-spectral methods, constructs an inter- 
polation of the magnetic field and then modifies the point values so 
that with the given interpolation scheme they produce a divergence- 
free continuous field. It is also possible to store and evolve point 
values of the magnetic vector potential, interpolate this vector po- 
tential, and find a value and derivatives of the magnetic field from 
this interpolation. This approach is used in the PENCIL CODlQ 
The vector potential approach always yields a magnetic field which 
is V • B = 0. Some of the disadvantages of this method are that it 
has the property that more than one vector potential configuration 
leads to the same magnetic field configuration, boundary conditions 
may be difficult to arrange, and compared to evolving B directly 
an extra level of spatial derivatives needs to be evaluated. 

The Smoothed Particle Method (SPH) attempts at MHD are 
notable in that SPH is a non-polynomial method used for approx- 
imating derivatives. In this context, in addition to the aforemen- 
tioned methods following the strategy of admitting V • B 7^ 0, 
a method based on the Euler angles formulation have been pro- 
posed by ,Price & BateH2007^ ; |Rosswog & PriceH2007^ w hich by 
construction yields V ■ B = , but Brandenburg 1 2010[ l has ob- 
served that this approach is not sufficient for realistic MHD as it 
severely constrains the allowed magnetic field geometries. Addi- 
tionally, ^rjce](|20T0| explored the use of the vector potential strat- 
egy in SPH, but found it to be unworkable. 
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In this paper, we describe a principle tiiat if adhered to allows 
point-value methods to evolve the magnetic field directly, while 
maintaining formally V ■ B = . Although throughout this paper 
we refer to magnetic fields, the principles and methods can apply 
to any vector field. 



2 A PRINCIPLE 

Since the discrete point values of the magnetic field do not have 
a defined derivative, the problem of V • B = lies entirely in 
the method used to produce the continuous representation of the 
magnetic field from which derivatives are taken. Thus, to pro- 
duce a V • i? = method, it is sufficient to define an interpola- 
tion (or quasi-interpolation) which is restricted to producing only 
V • B = fields. A concrete example of such a method is provided 
by divergence-free matrix valued radial basis function interpolation 
( |Narcowich & Ward|1994"||Low itzsch 2002| [2005] l. The following 
two sections of this paper are devoted to a summary of this tech- 
nique, and its use to construct finite-difference like operators. 

Radial Basis Function (RBF) Interpolation is an alternate 
method to polynomial basis interpolation for constructing functions 
which interpolate some discrete set of values. Instead of using a set 
of functions with a different polynomial form all mentored at the 
same place (a Taylor Series), RBF uses shifted versions of a one- 
parameter function. Further, these functions are shifted to be cen- 
tred on each interpolating point. For a set of scalar valued samples 
{xj, dj}f-i where Xj is the position of each point and dj is the 
value of the scalar field to be interpolated at that point, the RBF 
interpolant is of the form 

JV 

six) ^^i^iWx- Xj\\)Cj (I) 

where s{x) is the interpolant, i/; is a radial basis function, and Cj 
are a set of coefficients. These coefficients {cj}^i are such that 

s{xk) = dfe for alll ^ fc ^ iV. (2) 

Solving for these coefficients is done by solving the equation Gc — 
d where the matrix entries Gi,j = "l/dlxi — Xj\\). The remarkable 
ability of RBF interpolation is that if has certain properties, this 
equation has a unique solution for any set of points {xj} in any 
number of dimensions. The reader is encouraged to refer to | Wend-] 
|land ' (2005 1 and Buhmann (2003i for the mathematical details of 
the theory of radial basis function interpolation. 

Beyond scalar fields, it is possible to construct a RBF inter- 
polation for a vector field such as the magnetic field. If the RBF is 
chosen appropriately, this interpolation can be constrained to pro- 
duce V • B = 0. Given a set of point values {xj , dj where Xj 
is the position of each point and dj is the value of the vector field 
to be interpolated at that point, the interpolation is constructed in 
the form 

JV 

s(a;) = ^<E.(||a:-a;,||)c, (3) 

j=i 

where {cj}^i are such that 

8{xk) = dk for alll ^ fc sC iV (4) 
The matrix valued radial basis function <J> is constructed by 
^(a;) = {VV^ - V^/}V'(a;) (5) 
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Figure 1. The x vector field component of eq.^with e = 3.5, the y com- 
ponent is the same rotated 90 degrees. 

where ip is a scalar valued radial basis function and I is an identity 
matrix. If a numerical method is built using this representation for 
the magnetic fields, then the results will be free of V • B 7^ errors. 
This use of a V • B = interpolation basis is a general principle, 
it could apply to other classes of basis, and spectral basis. 



3 DEMONSTRATION 

In a manner similar to Taylor-series based finite difference stencils, 
we can build generalised finite difference stencils using radial basis 
function interpolation. The procedure is the same as for Taylor- 
series based finite differences, we interpolate the data at a local set 
of points, then take derivatives of the interpolant. Like with Taylor- 
series based finite differences, the resulting scheme will not in gen- 
eral conserve mass, linear momentum, or energy to machine pre- 
cision. These quantities will usually only be conserved to the level 
of the truncation error of the scheme. One can look to the body of 
work produced with the PENCIL CODE, a high order finite differ- 
ence scheme, to see examples of a successful approach based on a 
non-conservative method (ex:'Jo hansen et al^ ( |2007| l; |Babkovskaia,| 
Haugen & Brandenburg (201 1 1). Scalar value radial basis function 
finite difference stencils have been studied in |Bayona et al.|p010[ > 
for the case of the multiquadric radial basis function. The radial 
basis function finite difference approach (RBF-FD) has also been 
applied to convection-type PDEs in |Fomberg & Lehtol ( |20I0| l. 

To illustrate this construction, we must choose a radial basis 
function, in this case a Gaussian: 

■!/'(?■) — e'^^ (6) 

where e is a constant called the scaling factor. The scaling factor 
can be adjusted depending on the interpolation point distribution. 
Other radial basis functions can be used (see |Wendlandl|2005 1 or 
|Buhmann|2003| and recent results on the near equivalence of some 
common RBFs [Boyd 2010 1, but the Gaussian gives the simplest 
algebraic expressions in the following. 

To construct a divergence free matrix valued basis function 
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Figure 2. Upper. Alfven wave solutions at two resolutions, grid size N X N 
with 3x3 and 5x5 point stencils. Lower: Convergence of Alfven wave 
solution to analytical result. Results from 3x3 stencils plotted with + and 
from 5x5 stencils plotted with X . The solid line marks a slope of 1. The 
3x3 stencil error saturates at a larger value than the 5x5 stencil error. 



from tp{r), we apply equationjsjin two dimensions with = x'^ + 



(7) 



y , yielding: 

$11 = -(4e'y^ - 2e)e 
= 4e e ^ " ' 

$21 = 

$22 = — (4e X — 2e)e ^ " ' 



The combinations ("I'll, '&12) and ($21, '1>22) are divergence-free 
vector fields. Figure[T]shows the two components. One component 
resembles a dipole field in the x direction, and the other a dipole 
field in the y direction. In essence the method interpolates only 
V ■ B — fields because the field is built entirely from shifted 
and normalised versions of these dipole components. To build up 
an interpolation of some point-sampled field with these as the basis 
functions, it is necessary to solve the set of linear equations: 



For A'^ interpolation points the matrix A has entries: 

Al:JV,l:JV Ai,j = $ll(rij) 

An + 1:2N,1:N ^i.j = ^12{rij) 

Ai:N,N+l:2N — ^ Aij — $21 (^ij) 

^]V+1:2JV,]V+1:2JV — > Aij — ^22{rij) 



(9) 
(10) 
(11) 
(12) 



Each sub-matrix of A has entries corresponding to an entry in $. 
The vector d has entries: 



<il:JV 

<i]V+l:2JV 



di = Bi, 
di = Bi. 



Bxo 

ByO 



(13) 
(14) 



A is the interpolation matrix, and d is the values being interpo- 
lated. Bi^x and Bi.y are the components of the vector field being 
interpolated. Bxo and Byo are constant background field compo- 
nents, which may be chosen to be the field at the interpolation point 
where the derivatives are being calculated. This subtraction of the 
background constant component of the field increases the accuracy 
of the radial basis function approximation as this component is not 
in the space spanned by the interpolation basis. This background 
component is irrelevant to the V ■ B = constraint and to the 
determination of derivatives. The vector c is composed of the inter- 
polation coefficients in eq.[3] To find the derivatives of the interpo- 
lating function at point xo, we can simply evaluate the derivative 
of eq.|3] as 



ds 

dx 



N 

E 



a<E>| 

dx I 



(15) 



which yields the radial basis function estimate of the derivative 
at the point xq. This gives us a method of finding the derivatives 
of a V • B = magnetic field from point values. The interpola- 
tion points chosen can be arbitrary, but for the purposes of building 
finite-difference like derivative operators a set of nearest neighbour- 
ing points is natural choice. In the following, we demonstrate the 
use of 3 X 3 (9 point) and 5x5 (25 point) stencils, centred on xa, in 
two dimensions to solve the equations of magnetohydrodynamics. 

The equations solved are those for viscous, resistive, com- 
pressible isothermal magnetohydrodynamics in two dimensions: 



dp 
dt 

dpv 

at 



= -V ■ [pv) (16) 
= V ■V{pv)-VP + {V X B) X B + uV'^{pv) 

(18) 



with the equation of state P = c^p where p is the density, v is 
the velocity, P is pressure, B is the magnetic field, i' is the dy- 
namic viscosity, and 77 is the magnetic diffusivity. The equations are 
spatially discretized on an evenly spaced square grid with periodic 
boundary conditions. Spatial derivatives are estimated using the ra- 
dial basis function methods on 3 x 3 stencils as outlined above, 
and explicit time integration is performed with the forward Euler 
method. Both the derivatives of the scalar fields (p,P, components 
of v) and the vector field B are taken with scalar and vector RBF 
interpolations. This method is chosen so that the resulting code is 
as simple as possible to facilitate the reader's understanding. The 
source code in Python is available on the author's websit^ 

Since the method is V ■ B = by construction, we need only 



Ac = d 



(8) 
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Figure 3. Magnetised blast problem Left: Kinetic energy density Middle: Magnetic energy density Right: Mass density 



Table 1. Linear momentum errors in the magnetised blast problem 




CO 



11 1 ^' ' V' 1} ; I'lili ;1 



).00 



0.05 



0.10 
t 



O.IS 



Figure 4. Magnetised blast problem total mass error for two grid resolutions 
Af = 48 and Af = 96 
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Figure 5. Magnetised blast problem momentum errors for two grid resolu- 
tions Af = 48 and Af = 96 
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demonstrate that the solution converges as V ■ i? 7^ errors cannot 
occur. A suitable test is the evolution of a damped Alfven wave for 
finite V and rj, the analytical solution for which is given in |Chan-| 
|drasekhar| ( |1961[ >, section 39. The experimental convergence of the 
numerical solution to the analytical result for an Alfven wave is 
shown in Figure|2]with e = 1/64, u = ^ = 0.001. Note that the 
error saturates in this test for the 3x3 stencil. As the RBF interpo- 
lation used does not reproduce a first-order polynomial exactly, the 
approximation effectively stops improving below some grid spac- 
ing. To obtain further convergence, a larger stencil must be used. 
The Li error when the 5x5 point stencil is used decays at a rate of 
1.0, which is the limit set by the first order time integration scheme. 

As a further demonstration, the result of a magnetised blast 
problem, starting for an initial over density in the centre of the box 
is shown in Figure [3] The initial condition of this problem is = 
77 = 0.005, c. = 0.4082, p = 99e-((-o.5)^+(y-o.5)Vo.i2^)^ + 
V = 0, Bx = cos(27r/21). By = sin(27r/21) in a area 1x1 with 
a 96 X 96 grid with periodic boundary conditions. The output is 
shown at t = 0.2. The 3x3 stencil was used with e = 1/16. A 
similar problem is shown in |Stone et al.| ( |2008] l. We note the the in- 
termediate shock can be seen in the solution on the axis of the blast 
along the magnetic field jFerriere, Mac Low & Zweibel|199l| . The 
time history of the absolute error in total mass in the magnetised 
blast problem is shown in Figure |4] for two resolutions. The error 
is at the limit of numerical precision for all resolutions. Evidently 
the scheme used here is in fact mass-conserving even though it was 
not explicitly constructed to be so. The error in linear momentum 
is limited by truncation error, Figure|5]shows the time evolution of 
the absolute momentum error for two resolutions. Table [T] lists the 
momentum errors at t = 0.2 in the blast problem for four grid res- 
olutions A*'. The momentum error can be seen to converge towards 
zero as the resolution is increased. Energy conservation errors are 
not treated as the discussion is limited to isothermal magnetohy- 
drodynamics. 

The stability of this scheme is not explored here as the method 
is presented only as a demonstration of the use of these radial ba- 
sis function based derivative operators. The computational cost of 
a RBF based finite difference scheme for a derivative is the same as 
that for a polynomial based scheme with the same number of stencil 
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points, as the only difference is in tlie stencil coefficients. However, 
with RBF based schemes the well-motivated stencils may not be 
the same shape and size as for polynomial based finite differences 
- for example the square 3x3 stencil used here is not a popu- 
lar choice when used with polynomial based schemes. In contrast 
to the most directly comparable polynomial based finite difference 
scheme yielding V • i? = ( the vector potential method), the 
RBF based approach has the advantage of requiring fewer deriva- 
tive stencils to be computed as the magnetic field is obtained di- 
rectly not computed from derivatives of the vector potential. 



4 EXTENSIONS 

The method for constructing radial basis function based derivative 
approximations has no dependence on regularly placed points or 
the existence of mesh edges. Hence, these methods are easily used 
in a mesh-free context. Also, though in radial basis function inter- 
polation theory the interpolation points are chosen to be the radial 
basis function centres, in practise the approximation matrix can still 
be inverted if the interpolation points do not coincide with the ra- 
dial basis function centres. Furthermore, fewer radial basis func- 
tions can be used than interpolation points are specified - in this 
case a V • B = least squares approximation can be computed. 
The divergence-free basis used as an example in this work is not or- 
thonormal. If an orthonormal basis were specified, it would be pos- 
sible to preform divergence-free pseudo-spectral simulations with 
B, which may be of particular use in general relativistic MHD and 
force-free electrodynamics simulations. 

Other constraints beyond divergence-free can be placed on 
the vector field. For example, |LowitzscTi| |2005| observed that 
V X i? = type vector fields can be interpolated in a similar 
manner to shown here. This suggests the possibility to satisfy more 
complicated, though homogeneous, constraints. 



5 CONCLUSION 

In a point-value method, V ■ B = can be satisfied by the correct 
choice of interpolation scheme. Matrix-valued radial basis func- 
tions provide such an interpolations scheme. Finite-difference-like 
V ■ B = derivative operators can be constructed from matrix- 
valued radial basis function interpolations, and their use in the so- 
lution of magnetohydrodynamics problems has been demonstrated. 
Further exploration of the stability properties, accuracy, and com- 
putational cost of schemes based on these operators is warranted. 
The underlying principle of the choice of a V • B = interpola- 
tion also applies to pseudo spectral methods, and in general can be 
applied to any vector field where such a constraint is required. 
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